Efficient implementation of Gauss collocation and 
Hamiltonian Boundary Value Methods 

Luigi Brugnano* Gianluca Frasca Caccia^ Felice lavernaro''- 

, Devoted to celebrate the SO*'* birthday of John Butcher 

o 

Ph Abstract 

"^ In this paper we define an efficient implementation for the family of low-rank energy- 

C^ conserving Runge-Kutta methods named Hamiltonian Boundary Value Methods (HB- 

^ VMs), recently defined in the last years. The proposed implementation relies on the 

-^ particular structure of the Butcher matrix defining such methods, for which we can de- 

^^ rive an efficient splitting procedure. The very same procedure turns out to be automat- 

ically suited for the efficient implementation of Gauss-Legendre collocation methods, 
since these methods are a special instance of HBVMs. The linear convergence analysis 
of the splitting procedure exhibits excellent properties, which are confirmed by a few 
numerical tests. 
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2 1 Introduction 

• ^ The efficient numerical solution of implicit Runge-Kutta methods has been the subject of 

rS many investigations in the last decades, starting from the seminal papers of Butcher [23, 24] 

Cd (see also [25]). This aspect is even more relevant when dealing with geometric Runge-Kutta 

methods, that is, methods used in the framework of Geometric Integration where, usually, 

the discrete problems generated by the methods need to be solved to within full machine 

accuracy, in order not to waste the specific properties of the methods. 

In more details, in this paper we shall deal with the numerical solution of Hamiltonian 
problems, namely problems in the form, 

y' = JVH{y), y{to)=yoeR^^, (1) 
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where 

H{y) is the (scalar) Hamiltonian function defining the problem, and /„ the identity matrix 
of dimension m} Due to the skew-symmetry of matrix J one has 

l^H{y{t)) = VH{y{t)fy'{t) = VH{y{t)f JVH{y{t)) = 0, 

so that 

H{y{t)) = H{yo), t>h. 

For isolated mechanical systems, the Hamiltonian has the physical meaning of the total energy 
of the system, so that often the Hamiltonian is referred to as the energy. Its conservation 
is, therefore, a significant feature for the discrete dynamical system induced by a numerical 
method for solving (1): methods having this property are usually called energy- conserving 
methods. Among such methods, we are interested in the class of energy- conserving methods 
named Hamiltonian Boundary Value Methods (HBVMs) [9, 10, 11, 8, 12, 13, 14, 15, 16, 4, 5, 
6, 3], which have been recently devised starting from the idea of discrete line integrals, defined 
in [29, 30, 32] (see also [31]). For such methods, the discrete problem can be conveniently 
posed in a suitable form which can be exploited to derive efficient implementation strategies, 
as was done in [12]. Here we further improve on such results, by proposing and analysing 
an iterative procedure based on the particular structure of the discrete problem. As a by 
product, an efficient implementation of Gauss-Legendre collocation methods is also obtained. 
Indeed, these latter methods may be interpreted as a particular instance of HBVMs. The 
proposed procedure is strictly related to that recently devised in [7] for Radau HA collocation 
formulae, though the two approaches are substantially different. 

With this premise, the paper is organized as follows: in Section 2 we describe the structure 
of the discrete problem generated by HBVMs, along with the way of solving it, as done so far; 
in Section 3 we introduce the new iterative procedure, which is based on a suitable splitting; 
in Section 4 we study the convergence properties of the new iteration, also comparing it with 
known existing ones; in Section 5 a few numerical tests are reported; at last, a few conclusions 
are contained in Section 6. 

2 Discrete problem induced by HBVMs 

We now recall the basic facts about HBVMs, and derive the most efficient formulation of the 
generated discrete problems. Let us assume, for sake of brevity, to = in (1), and consider the 
approximation of the problem over the interval [0, /i], which will clearly concern the very first 
application of a given numerical method. Let us then consider the orthonormal polynomial 
basis, on the interval [0, 1], provided by the shifted and scaled Legendre polynomials {Pj}'. 

degPj = i, / Pi{x)Pj{x)dx = 5ij, ^i,j>0, (3) 



^In the following, when the size of the identity is not specified, it can be deduced from the context. 



where 5ij is the usual Kronecker symboL Under suitable mild assumptions on the Hamiltonian 
function H ^ the right-hand side of the differential equation (1) can be expanded along the 
considered basis, thus giving 

y\ch) = Y,lAy)PAc), cG [0,1], (4) 

i>o 

where 

7,(1/) = / JV//(y(r/i))P,(r)dr, j > 0. (5) 

Jo 

By imposing the initial condition, the solution of this problem is formally obtained by 

y{ch) = yo + hJ2lAy) [ PA^)dx, c E [0, 1]. (6) 

In order to derive a polynomial approximation a of degree s to (6), we consider the following 
approximated ODE-IVPs: 

s-l 

a'(c/i) = 5^7,(a)P,(c), CG [0,1], a(0)=i/o, (7) 

j=0 

where 7j(o-) is defined according to (5), by formally replacing y by a. Consequently, the 
approximation to (6) will be given by 

a{ch) = yo + hy2lM) / Pjix)dx, c G [0, 1]. (8) 

,=0 Jo 

Assume now that the Hamiltonian function is a polynomial of degree u. Consequently, the 
(unknown) vector coefficients {7,(0")} can be exactly obtained by using a quadrature formula 
defined at the Gaussian abscissae {ci, . . . , Ck}, i.e., 

Pk{c,) = 0, i = l,...,k, (9) 

and corresponding weights {61, . . . , bk},"^ 

k 

7,M = Yl hJVH{a{cM)Pj{c^), J = 0, . . . , s - 1, (10) 

provided that 

.<f. (U) 

By setting 

Yi = a{c^h), i = l,...,k, (12) 



^Hereafter, we shall always assume this choice. 



and considering that the new approximation is given by 



yi 



/■I *; 

a{h) = yo + h JVH{a{Th))dT = y^ + hY, hiJVH{Y,), 
Jo i=l 



one then obtains the following A;-stage Runge-Kutta method, 

A 



(13) 



where, as usual, b,c ^ 'R'' are the vectors containing the weights and the abscissae, respec- 
tively, and (see, e.g. [12, 13, 14]) 



with 



A = Vs+iXsV^n G 



/ Po(ci) ... P.-i(ci) \ 



•tkxk 



(14) 



Vr 



( \ -6 

6 



nfcxr 



r = s, s + 1 



1" \ -^1 



(15) 



X 



^r-l(Cfc) / 



is-X 



V 



X 



o...oe. 



ps+lxs 



(16) 



(17) 
fl8) 



is I 

6 = (2V422-l)"\ ^ = l,...,s, 
^ = diag(6) e M^^^ 

According to [11], we give the following definition (see also [9]). 

Definition 1 (13)-(18) define a HBVM{k,s) method. 

The following properties [11, 14] elucidate the role of two indices k (number of ascissae) and 
s (degree of the underlying polynomial a) in the previous definition. 

Theorem 1 For all k > s, a HBVM{k, s) method: 

• has order 2s, that is: 

y, - y(h) = Oih'^^'y, 

• is energy conserving for all polynomial Hamiltonians of degree v satisfying (11); 

• for general non-polynomial (but suitably regular) Hamiltonians, one has: 

H{y,)-Hiyo)=Oih''^'). (19) 



Remark 1 From (19) one deduces that a HBVM{k, s) method is practically energy- conserving 
also in the case of non-polynomial Hamiltonians, provided that k is large enough. Indeed, it 
is enough to approximate the involved integrals to within round-off' errors. 

For sake of completeness, and for later reference, we also prove the following result, which 
actually shows that HBVM(A;, s) methods can be regarded as a generalization of the s-stage 
Gauss-Legendre collocation formulae [11]. 

Theorem 2 HBVM{s,s) coincides with the s-stage Gauss-Legendre collocation method. 

Proof. When k = s, from (9) and (14)-(18) one obtains: 

p,+i = {Vs 0) ^ A = PsXsVjn. 

Moreover, 

vjnvs = Is, 

since, by setting Cj the ith unit vector, 

k „1 

ej {VjnVs) e, =J2beP^-l{ce)PJ-lice) = / Pi-i{x)Pj_i{x)dx = 6^j, 
due to the fact that the quadrature is exact, and because of (3). Consequently, 

vjn = v;' 

and, then, 

A = VsXsV;\ (20) 

which is the W^-transformation [26] of the s-stage Gauss-Legendre collocation method. D 

If we set y the (block) vector with the internal stages (12) and e = (1, . . . , 1)-^ G M'^, the 
discrete problem generated by a HBVM(/c, s) method is given by 

y = e®yo + hA®JVH{y), (21) 

which is a nonlinear system of (block) dimension k.^ However, in view of (19), k is likely to be 
much higher than s and, consequently, such a formulation is in general not recommendable. 
To derive a more efficient formulation, let us set 7 the (block) vector containing the 
coefficients defining the polynomial a in (8), thus obtaining: 

7 = Vjn ® JVH{y), y = e®yo + hTs+iX,, ® /7- 

Combined together, such equations provide us with the following discrete problem, 

F{j) = J- Vjn ^.JVH(e^yo + hTs+iX^ 17) = 0, (22) 

■'Here VH{y) is the block vector whose entries are given by the gradient of H evaluated at the k stages. 
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whose (block) size is always s, independently of k. One easily checks that the simplified 
Newton iteration, applied for solving (22) consists in the following iteration [12]: 

solve: [l-hXs®JV^Ho]A^ = -F{Y) (23) 

y+^ = y+A^ £=o,i,..., 

where V^Hq is the Hessian of H{y) evaluated at yo. Consequently, the bulk of the computa- 
tional cost is due to the factorization of the matrix 

having dimension 2sm x 2sm. In the next section, we shall see how to efficiently solve the 
iteration (23). 

3 The new splitting procedure 

The iteration (23) is similar in structure to the simplified-Newton iteration applied to the 
original system (21), for which a number of splitting procedures have been devised (see, 
e.g., [27, 28, 1, 2, 7, 17, 18, 21, 19, 20]). However, the triangular splitting iteration defined 
in [27, 28], and the modified triangular splitting iteration defined in [1] turn out to be not 
effective for (23), due to the particular structure of the matrix Xs (see (16)). Conversely, the 
blended iteration defined in [2, 17, 18, 21, 19, 20] turns out more appropriate, as is shown 
in [12]. We here shall devise a different iterative procedure, which appears to be even more 
favourable. This is the subject of the remaining part of this section. The main idea is similar 
to that devised in [7] for Radau HA collocation methods, even though the framework and 
the overall details (and results) are different. 
To begin with, let us consider the polynomial 

s-l 

7(c) = 5^P,-(c)7,(a), ceR, (24) 

j=0 

and introduce a new set of (block) unknowns, 

s-l 

% = ^Pj{ci)-fj{a), i = l,...,s, (25) 

j=0 

defined as the evaluation of (24) at the set of auxiliary abscissae 

ci < • • • < Cs- (26) 



Introducing the (block) vector 



7 



/7i\ 



V7. / 



(27) 



and the matrix 

V = ( P,_i(q) ) G M^x^ (28) 

we can recast (25) in matrix notation as 

j = V^Ij. (29) 

In terms of the new unknown vector 7, the simphfied Newton iteration (23) reads: 

solve : MoA^ = -V ® I FiV'^ ® I Y) = 4 (30) 

f+i = f + M £ = 0,1,..., 

where 

MQ = I-h {vX,V-^\ (g) JV^Ho = I-hA(g) JV^Hq. (31) 

The key idea is that of choosing the abscissae (26) such that the Grout factorization of A 

A = LU, (32) 

has the diagonal entries of L all equal. In such a case, by following the approach of van der 
Houwen et al. [27, 28], the iteration (30) is replaced by the inner-outer iteration 



solve : 



I ~hL® JV^Ho 



r = 0,l,...,z/-l, (33) 



7^+^ = 7^ + A^'^ £ = 0,1,.... 

In particular, since A^''^ = 0, the choice z/ = 1 corresponds to the approach used by van 
der Houwen et al. to devise PTIRK methods [27], whereas, if u is large enough to have full 
convergence of the inner-iteration (the one on r), then the outer iteration is equivalent to 
(30). Clearly, all the intermediate possibilities can be suitably considered. 

That the diagonal entries of the factor L are all equal to a given value, say ds, has the 
obvious advantage that one only needs to factor the matrix 

/ - hdsJV^Ho e M2mx2m^ 

Concerning dg, the following result holds true. 

Theorem 3 Assume that the factorization (32) is defined and that the factor L has all its 
diagonal entries equal to dg. Then, with reference to (17), one has: 



'n,--i Cl-i ) if s is even, 

ds={ \ (34) 

' nil! a , ^fsis odd. 



s 



Proof. Assume that (31)-(32) hold true. Then 

det(X,) = deiiVXsV^^) = det(i) = det{LU) = det(L) = d^, 
since U has unit diagonal and all the entries of L are equal to dg- Consequently, 

ds = Vdet(X,). 
The thesis then follows by considering that, from (16), 

det(Xi) = ^, det(X2) = ei, 
and, by and applying the Laplace formula, one obtains: 

■ n»-i ^2i-i 5 if s is even, 

det(X,) = <( '''-]:'' (35) 

^miki., if sis odd. 



By virtue of the previous result, in order to compute the auxiliary abscissae (26), we have 
symbolically solved the following set of equations, which is obviously equivalent to requiring 
that the factor L has the diagonal entries equal to each other: 

det(i^+i) =4det(i£), i=l,...,s-l. (36) 

Here Af denotes the principal leading submatrix of order i of A, and dg is given by (34). 

Remark 2 We observe that the auxiliary abscissae (26) are s, whereas the conditions (36) 
are s — 1. This means that a further condition can be imposed on the abscissae. We have 
chosen it in order to (approximately) minimize the maximum amplification factor of the 
iteration, as described in Section 4- 

The obtained results are listed in Table 1, for s = 2, . . . , 6, from which one sees that in all 
cases the abscissae are distinct and inside the interval [0, 1]. 

We emphasize that, for any given s, the distribution of the abscissae q is independent of 
k and so is the factorization (32) of the matrix A whose computation is responsible of the 
bulk of the computational effort during the integration process. This property has a relevant 
consequence during the implementation phase of this class of methods. In fact, one can 
conjecture a procedure to advance the time that dynamically selects the most appropriate 
value of k. Depending on the specific problem at hand and the configuration of the system 
at the given time, one can easily switch from a symplectic to an energy preserving method 
by choosing k = s (Gauss method) or k > s, respectively. 

4 Convergence analysis and comparisons 

In this section we briefly analyze the splitting procedure (33) according to the linear analysis 
of convergence in [27] (see also [20]). In such a case, problem (1) becomes the celebrated test 
equation 

y' = Ay, y{to) = Vo- (37) 



Table 1: Auxiliary abscissae (26) for the HBVM(A;, s) and s-stage Gauss method, s = 2, 
and the diagonal entry ds (see (34)) of the corresponding factor L. 



.,6, 



s = 2 


Cl 
C2 


0.26036297108184508789101036587842555 
1 


d2 


0.28867513459481288225457439025097873 



Cl 0.15636399930006671060146617869938122 
C2 0.45431868644630821020177903150137523 
C3 0.948 



4 0.20274006651911333949661483325792675 



Cl 0.10980739789315927030006609946548674 

C2 0.33170147176242312021431664090897753 

C3 0.34895222925320098805117573706272031 
C4 0.7186 



di 0.15619699684601279005430416526875577 



Cl 0.084665152651785050300050089866691763 

C2 0.257176246353319497655094082620443041 

C3 0.406997899472141237569051120583722535 

C4 0.543482236497480925424020865878224100 
C5 0.8858 



4 0.12702337351164258963093490787943281 



6 



Cl 0.067421939209393759076155059271608432 

C2 0.199433361269530107247831326370752698 

C3 0.375498823199182884162061752054007488 

C4 0.44562087499073615218681054513866021 

C5 0.65709214779914978056665114753638658 

Cg 1 



4 0.10702845478806509529222890981996019 



By setting, as usual, q = hX, one then obtains that the error equation associated with the 
outer iteration (33) is given by 

e,+i = Z{qyee, Z{q) := q{I - qL)~'L{U - /), £ = 0, 1, ... , (38) 

where e^ = A^'" — A is the error vector at step i and Z{q) is the iteration matrix induced by 
the sphtting procedure. This latter will converge if and only if its spectral radius, 

p(g):=p(Z(g)), 

is less than 1. The region of convergence of the iteration is then defined as 

B = {qeC : p(g) < 1}. 

The iteration is said to be A-convergent if C^ C D. If, in addition, the stijf amplification 
factor, 

p~ := lim p{q), 

(jr— >oo 

is null, then the iteration is said to be L-convergent. Clearly, A-convergent iterations are 
appropriate for A-stable methods, and L-convergent iterations are appropriate for L-stable 
methods. In our case, since 

Z{q)^{U-I), q^oo, (39) 

which is a nilpotent matrix of index s, the iteration is L-convergent if and only if it is A- 
convergent. Since the iteration is well defined for all q G C~ (due to the fact that the diagonal 
entry of L, dg, is positive, as was shown in (34)) and p(0) = 0, A-convergence, in turn, is 
equivalent to require that the maximum amplification factor, 

P* = maxp(ix), 

is not larger than 1. Another useful parameter is the nonstiff amplification factor, 

p = p(L(f/-/)), (40) 

that governs the convergence for small values of q since 

p(g) ~ p\q\, for q ^ 0. 

Clearly, the smaller p* and p, the better the convergence properties of the iteration. 

Remark 3 As anticipated in Remark 2, the additional condition imposed on the auxiliary 
abscissae (26) is that of (approximately) minimizing the maximum amplification factor p* of 
the iteration, while fulfilling the conditions (36). 

In Table 2 we list the nonstiff amplification factors and the maximum amplification factors 
for the following L-convergent iterations applied to the s-stage Gauss-Legendre methods: 

(i) the iteration obtained by the original triangular splitting in [27]; 
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Table 2: Amplification factors for the triangular splitting in [27], the modified triangular 
splitting in [1], the blended iteration in [17], and the splitting (33), for the s-stage Gauss- 
Legendre formulae. The last two cases coincide with those for the HBVM(A;, s) methods, 
k> s. 





(i): triangular 


(ii): triangular 


(iii): 


blended 


(iv): triangular 




sphttin} 


y in [27 


splitting in [1] 


iteration in [17 


splitting (30) 


s 


P 


P* 


P P* 


P 


P* 


P P* 


2 


0.1429 


0.0833 


0.1340 0.0774 


0.1340 


0.0774 


0.1340 0.0774 


3 


0.3032 


0.1098 


0.2537 0.0856 


0.2765 


0.1088 


0.2536 0.0870 


4 


0.4351 


0.1126 


0.3492 0.0803 


0.3793 


0.1119 


0.2656 0.0696 


5 


0.5457 


0.1058 


0.4223 0.0730 


0.4544 


0.1066 


0.3275 0.0534 


6 


0.6432 


0.0973 


0.4861 0.0702 


0.5114 


0.0993 


0.4243 0.0878 



(ii) the iteration obtained by the modified triangular splitting in [1]; 

(iii) the blended iteration obtained by the blended implementation of the methods, as defined 
in [17]; 

(iv) the iteration defined by (33). 

We recall that the scheme (i) (first column) requires s real factorizations per iteration, whereas 
(ii)-(iv) only need one factorization per iteration. From the parameters listed in the table, one 
concludes that the proposed splitting procedure is the most effective among all the considered 
ones. 

Remark 4 For sake of accuracy, we stress that, when dealing with the actual implementation 
of HBVM{k,s) methods, only the blended iteration [12] and the newly proposed one (33) can 
be considered, whereas the triangular splitting defined in [27] and its modified version [1] turn 
out to be not effective, as was pointed out at the beginning of Section 3. Consequently, in 
such a case, one has to consider only the last two group of columns in Table 2. 



5 Numerical Tests 

In this section, we report a numerical example taken from [12], namely the Hamiltonian 
problem defined by the Hamiltonian 

H{q,p) =p' + /32g2 + ^2^q + ^fn ^4^) 

with the following parameters (slightly modified, with respect to those considered in [12]): 

/3 = 12, a = 4, n = A. (42) 



11 



Table 3: Number of fixed-point, blended iterations, and splitting iterations (- means no 
convergence) required for each level curve of problem (41)-(43), by using the HBVM(8,2) 
and HBVM(12,3) methods with stepsize h = 10^^, along with the maximum error in the 
numerical Hamiltonian. 





HBVM(8,2) 




HBVM( 


12,3) 






Fixed-point 


Blended 


Splitting 




Fixed-point 


Blended 


Splitting 




i 


iterations 


iterations 


iteration 


errH 


iterations 


iterations 


iteration 


errH 





9191 


6957 


4814 


2.9 


10-^5 


8544 


7398 


4681 


3.1 


10-15 


1 


16835 


12151 


8648 


1.6 


10-14 


14684 


12486 


8086 


7.4 


10-15 


2 


23907 


16344 


12293 


6.6 


10-14 


19004 


16860 


11095 


3.2 


10-14 


3 


33396 


21854 


17089 


3.2 


10-14 


23854 


22554 


15415 


9.8 


10-14 


4 


50163 


30762 


23744 


5.4 


10-14 


31263 


31463 


22445 


9.3 


10-14 


5 


- 


48109 


35994 


6.5 


10-14 


42598 


45747 


34885 


1.5 


10-13 


6 


- 


- 


51961 


9.9 


10-13 


- 


71105 


49239 


2.9 


10-13 


7 


- 


- 


- 


- 


- 


- 


66500 


2.3 


10-13 



Consequently, the Hamiltonian is a polynomial of degree z/ = 8 and, from (11), one has 
that HBVM(8,2) and HBVM(12,3) methods are energy-conserving, with orders 4 and 6, 
respectively. We integrate the orbits starting at the initial points 



(go, po ) = ( ^ + 0.1, -(^ + 0.1) ) 



0, 



7. 



(43) 



These orbits are closed (see Figure 1) and in [12] it has been shown that, for their numeri- 
cal integration, energy-conserving methods are more appropriate than symplectic integrators 
(like Gauss-Legendre collocation methods). We use a constant stepsize h = 10-3 for inte- 
grating the eight trajectories for 103 steps. In Table 3 we list the results obtained by the two 
methods, in terms of required iterations for solving the generated discrete problems, when 
using: 

• a fixed-point iteration; 

• a blended iteration; 

• the splitting iteration (33) with 2 inner iterations. 

The choice of 2 inner iterations in (33) makes the cost of one outer iteration comparable to 
that of one blended iteration, provided that (33) is implemented as suggested in [7].'* 

From Table 3 we see that the proposed splitting iteration (33) appears to be the most 
reliable among the considered ones and, in particular, it is competitive with the blended 
implementation of HBVMs. Similar results are obtained when considering Gauss-Legendre 
collocation methods, even though the qualitative behavior of the computed solution may be 
not correct when moving away from the origin (see [7]). 



*The implementation details of the new procedure will be the object of a future work. 



12 



10 



4 

2 

Q. 

-2 
-4 



-10 




Figure 1: Level curves for problem (41)-(43). 

6 Conclusions 

In this paper we have defined an efficient iterative procedure for solving the discrete problems 
generated by the application of HBVM(A;, s) methods, a class of energy-conserving methods 
for polynomial Hamiltonian dynamical systems. The proposed implementation turns out to 
improve over that proposed in [12], as is confirmed by a few numerical tests. In particular, it 
also applies to Gauss-Legendre formulae and the resulting linear convergence analysis shows 
that the proposed iterative procedure is the most effective, among those based on suitable 
splittings of the corresponding Butcher array of the methods. 
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